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ABSTRACT 

Context. The study of pre-stellar cores (PSCs) suffers from a lack of undepleted species to trace the gas physical properties in their very dense 
inner parts. 

Aims. We want to carry out detailed modelling of NiH + and N2D + cuts across the LI 83 main core to evaluate the depletion of these species 
and their usefulness as a probe of physical conditions in PSCs. 

Methods. We have developed a non-LTE (NLTE) Monte-Carlo code treating the ID radiative transfer of both N2H + and N2D + , making use 
of recently published collisional coefficients with He between individual hyperfine levels. The code includes line overlap between hyperfine 
transitions. An extensive set of core models is calculated and compared with observations. Special attention is paid to the issue of source 
coupling to the antenna beam. 

Results. The best fitting models indicate that i) gas in the core center is very cold (7± 1 K) and thermalized with dust, ii) depletion of N^H + 
does occur, starting at densities 5-7 x 10 5 cnr 3 and reaching a factor of 6+" in abundance, iii) deuterium fractionation reaches ~70% at the 
core center, and iv) the density profile is proportional to r _I out to ~4000 AU, and to r 2 beyond. 

Conclusions. Our NLTE code could be used to (re-)interpret recent and upcoming observations of N2H + and N2D + in many pre-stellar cores 
of interest, to obtain better temperature and abundance profiles. 

Key words. ISM: abundances - ISM: molecules - Radiative transfer - ISM: structure - ISM: individual: L183 - Line: formation 



1. Introduction 

Understanding star formation is critically dependent upon the 
characterisation of the initial conditions of gravitational col- 
lapse, which remain poorly known. It is therefore of prime im- 
portance to study the properties of pre-collapse objects, the so- 
called pre-stellar cores (hereafter PSCs). As bolometers and 
infrared extinction maps have been unveiling PSCs through 
their dust component, it has become clear that depletion of 
molecules onto ice mantles is taking place inside these cores, 
preventing their study with the usual spectroscopic tools. In 
most PSCs, very few observable species seem to survive in the 
gas phase in the dense and cold inner parts, namely N2HP, NH3, 
H2D + and their isotopologues (e.g. Tafalla et al. 2002). In a few 
cases, it is advocated that even N-bearing species also deplete 
(e.g. B68: Bergin et al. 120021 L1544: Walmsley et al. |2004l 
L183: Pagani et al. |2005l hereafter PPABC). 



Send offprint requests to: L. Pagani 

* Based on observations made with the IRAM 30-m and the CSO 
10-m. IRAM is supported by INSU/CNRS (France), MPG (Germany), 
and IGN (Spain). 



Among the above three species, H2D + is thus the only 
one not to deplete. However, it has only one (ortho) transi- 
tion observable from the ground, that moreover requires ex- 
cellent sky conditions. The para form ground transition at 1.4 
THz should not be detectable in emission in cores with Tkin < 
10 K. Therefore H2D + is useful for chemical and dynamical 
studies, but brings little information on gas physical conditions. 
NH3 inversion lines at 23 GHz can provide kinetic temperature 
measurements as long as higher lying non-metastable levels 
are not significantly populated (Walsmley & Ungerechts 1983 ). 
However, because the critical density of the NH3 (1,1) inver- 
sion line is only ~2000 cm~ 3 , this tracer may have substantial 
contribution from external, warmer layers, not representative 
of the densest parts of PSCs. The third species, N2H + , appears 
very promising : it has the strong advantage of having mm tran- 
sitions with critical densities in the range 0.5-70 x 10 6 cm -3 
and intense hyperfine components for the (J: 1-0) line in typi- 
cal PSC conditions. The ratio of hyperfine components gives an 
estimate of the opacity and excitation temperature of the line. 
Fitting both the hyperfine ratios and the (J: 1-0) to (J:3-2) ra- 
tio with a rigorous NLTE model should thus bring strong con- 
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Fig. 1. N2H + observations (grey lines) compared to our best radiative transfer model (see FigH^-b). Each column corresponds to 
a spatial offset (indicated above the top box) from the main dust peak. The top row contains the CSO 10-m data, the other rows 
the IRAM 30-m data. The (J: 1-0) line is split between the bottom row (central triplet) and the next row ('red' triplet); the isolated 
'blue' component is not shown. The rightmost column represents the fit for the central position with a NLTE code without overlap 
for the same cloud parameters. 



straints on both Tk; n and n(H2). Collisional coefficients (with 
He) between individual hyperfine sublevels have become avail- 
able recently (Daniel et al. 2005). However, current excitation 
models (Daniel et al. 2006a) do not take into account line over- 
lap, which limits their accuracy. Another question that remains 
is whether N2H + is indeed able to probe the central core re- 
gions, despite the depletion effects which have been reported. 

In this paper, we introduce a new NLTE Monte-Carlo ID 
codf] treating N2H + and N2D + radiative transfer with line over- 
lap, and apply it to detailed analysis of the main PSC in LI 83, 
a clear-cut case of N2H + depletion (cf. PPABC). We demon- 
strate the capability of the model to constrain physical condi- 
tions inside the PSC (temperature, density profile, abundance 
and depletion, deuterium fractionation). In particular we show 
that the gas is very cold at the core center and thermalized with 
the dust, and that N2D + appears a very useful tracer of physical 
conditions in the innermost core regions. 

2. Observations and analysis 



The main dust core in L183 (Pagani et al. 2004) was observed 
at the IRAM 30-m telescope in November 2003, July 2004 and 
October 2006. Spectra were taken on a 12" grid in an East- 
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West strip across the core (centered at (2(2000) = 15h54m08.5s, 
5(2000)= -2°52'48"). Symmetric eastern and western spectra 
were averaged out to ± 48" to give a more representative radial 
profile. A small anti-symmetric velocity shift of a few tens of 
m s _1 was noticed on either sides of the PSC center at distances 
beyond ± 30", possibly indicative of rotation (see Section [3~2l . 
This shift was compensated for when averaging eastern and 
western spectra, in order not to artificially broaden the lines. 
Beyond 48", only eastern positions are considered, as the west- 
ern side is contaminated by a separate core (Peak 3; Pagani et 
al. 2004). Lines of N 2 H + (J: 1-0), (J:3-2) and of N 2 D + (J:l- 
0), (J:2-l) and (J:3-2) were observed in Frequency switching 
mode, with velocity sampling 30-50 ms _I and T sys ranging 
from 100 K at 3mm up to 1000 K at 1mm. The N 2 H + (J:2-l) 
line at 186 GHz was not observed, as it lies only 3 GHz away 
from the telluric water line at 183 GHz, hence its usefulness 
to constrain excitation conditions would be limited by calibra- 
tion sensitivity to even tiny sky fluctuations. The problem does 
not apply to N2D + which lies a factor of 1 .2 lower in frequency. 
Spatial resolution ranges from 33" at 77 GHz to 9" at 279 GHz. 
Additional CSO 10-m observations of N 2 H + (J:3-2) were ob- 
tained in June 2004 at selected positions. Observations were 
done in Position Switch mode with the high resolution AOS 
(48 kHz sampling) and a T sys around 600 K. A major interest 
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Table 1. Spectral parameters of each line observed towards the PSC center position (ff(2000) = 15h54m08.5s, 5(2000)= 
-2°52'48") with V L sr = 2.3672(2) kms -1 (from the NH3 (1,1) measurement): noise level, reference observing frequency, 
intrinsic linewidth of hyperfine transitions and total opacity (both from the CLASS HFS routine), relative velocity and intensity 
of the main hyperfine groups (from gaussian fits). Values in parentheses are lcr uncertainties on the last digit. For comparison, 
the ratios of hyperfine groups given by their statistical weights are listed in italics. The detailed hyperfine structure of the 3 
species is described in Caselli et al. ( fl995T l for N 2 H + , Dore et al. d20Q4b for N 2 D + and Kukolich JT967b for NH 3 



Noise (lcr) : 
Ref . frequency : 
Intrinsic FWHM 
Total opacity : 

Hyperfine 
Group 



N 2 H + (J : 1 - 0) 

28 mK 
93.173764 GHz 
0.195(1) kms" 1 

21.8(5) 



VrA. 

(kms- 1 ) 



T R 

(K) 



N 2 H + (J : 3 - 2)(IRAM) 
59 mK 
279.511832 GHz 
0.18(2) kms -1 
10(2) 



Area T R Ratio V rc i 

(Kkms- 1 ) obs weights (kms -1 ) 



1 R 

(K) 



Area 
(Kkms- 1 ) 



N 2 H + (J : 3 - 2)(CSO) 
64 mK 
279.511832 GHz 
a 0.28(4)kms-' 
2(1) 



Kel. 

(km s _1 ) 



T R Area 
(K) (Kkms -1 ) 



1 


-8.009(7) 2.05(3) 


0.47(3) 


0.932 


0.428 


-2.58(1) 0.17(6) 0.013(5) 


-2.63(2) 0.20(6) 0.025(8) 


2 


-0.611(8) 1.81(3) 


0.42(3) 


0.822 


0.428 


0.00(1) 0.50(6) 0.29(2) 


0.00(1) 0.67(6) 0.37(2) 


3 


0.000(7) 2.20(3) 


0.60(3) 


1 


1 


4.70(8) 2.25(6) 0.014(5) 


4.68(3) 2.12(6) 0.028(9) 


4 


0.956(7) 2.08(3) 


0.54(3) 


0.945 


0.714 






5 


5.546(7) 2.12(3) 


0.47(3) 


0.964 


0.428 






6 


5.983(8) 2.04(3) 


0.53(3) 


0.927 


0.714 






7 


6.94(1) 1.26(3) 


0.24(3) 


0.573 


0.143 








N 2 D + (J : 1 - 


0) 






N 2 D + (J : 2 - l) b 


N 2 D + (J : 3 - 2) 


Noise (lcr) : 


31 mK 








26 mK 


43 mK 


Ref. frequency : 


77.109616 GHz 






154.217182 GHz 


231.321917 GHz 


Intrinsic FWHM : 


0.226 (3) kms- 1 






0.185 (2) kms- 1 


0.18 (2) kms-' 


Total opacity : 


4.7(3) 








4.9(1) 


1.5(7) 




Vra. T R 


Area 


Tr 


Ratio 


Frei. T R Area 


Vrei. T R Area 




(kms" 1 ) (K) 


(Kkms- 1 ) 


obs 


weights 


(kms- 1 ) (K) (Kkms- 1 ) 


(kms- 1 ) (K) (Kkms- 1 ) 



1 -9.697(8) 0.77(3) 

2 -0.763(8) 0.76(3) 

3 0.000(7) 1.10(3) 

4 1.146(8) 0.92(3) 

5 6.65(2) 0.55(3) 

6 7.19(1) 0.90(3) 

7 8.34(2) 0.34(3) 



0.15(1) 


0.7 


0.428 


0.15(1) 


0.691 


0.428 


0.33(2) 


1 


1 


0.25(2) 


0.836 


0.714 


0.16(2) 


0.5 


0.428 


0.23(2) 


0.818 


0.714 


0.05(2) 


0.309 


0.143 



-5.35(2) 0.16(3) 0.137(5) 

-0.759(8) 0.21(3) 0.080(3) 

0.000(1) 1.37(3) 0.513(3) 

0.646(4) 0.40(3) 0.139(3) 

2.5(1) 0.08(3) 0.02(2) 

2.97(2) 0.43(3) 0.08(2) 

3.57(5) 0.26(3) 0.16(3) 



-3.26(1) <0.12 

0.000(5) 0.72(4) 0.186(8) 

0.46(2) 0.19(4) 0.06(1) 

2.47(3) 0.11(4) 0.028(7) 



Noise (lcr) : 
Ref. frequency : 
Intrinsic FWHM 
Total opacity : 



NH 3 (1, 1) 
74 mK 
23.6944957 GHz 
0.195(1) kms-' 
24.2(4) 



Vrei. T R Area i R K.auu v rs [ 

(kms- 1 ) (K) (Kkms-') obs weights (kms-') 



T R Ratio 



NH 3 (2, 2) 
76 mK 
23.7226332 GHz 
0.20(1) kms- 1 
0.1(7) 



V„ 



Area 
(K) (Kkms- 1 ) 



T* 



1 -19.504(6) 2.62(7) 0.89(4) 0.922 0.363 -0.005(5) 0.73(8) 0.152(7) 

2 -7.814(5) 2.68(7) 0.69(3) 0.944 0.273 

3 -7.255(8) 1.95(7) 0.57(3) 0.687 0.182 

4 -0.17(1) 2.84(7) 1.15(9) 1 1 

5 0.30(2) 2.78(7) 1.0(1) 0.979 0.636 

6 7.465(9) 2.60(7) 0.72(5) 0.915 0.303 

7 7.89(1) 1.91(7) 0.48(5) 0.673 0.152 

8 19.318(9) 2.48(7) 0.62(4) 0.873 0.242 

9 19.85(1) 1.72(7) 0.42(5) 0.606 0.121 



a CSO observations are done with a resolution around 0.1 kms 1 (or more) which explains this larger value 

b The hyperfine structure is too complex to be detailed completely. Only the main groups of hyperfine transitions are given 
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to observe the (J:3— 2) line with the CSO is that the beam size 
and efficiency is very similar to the 30-m values for the (J:l- 
0) line, thus almost canceling out beam correction errors in the 
comparison between the two. It is thus a useful constraint for 
radiative transfer modelling, even though the higher resolution 
(9") 30-m (J:3-2) observations will be essential to constrain 
abundances and physical conditions in the innermost core re- 
gions. Standard reduction techniques were applied using the 
CLASS reduction package^. 

Complementary observations of NH3 (1,1) and (2,2) inver- 
sion lines towards the PSC center were also obtained at the 
new Green Bank 100-m telescope (GBT) in November 2006 
with velocity sampling of 20 m s" 1 and a typical T sys of 50 K, 
in Frequency Switch mode. The angular resolution (~35") is 
close to that of the 30-m in the low frequency (J: 1-0) N2D + line. 
The main beam efficiency of the GBT at 23 GHz is between 
0.95 and 1, hence no correction was applied to put the spectra 
on T^ scale. 

Table Q] summarizes the main observational parameters of 
the lines observed towards the core center: noise level of 
the spectrum (rms), rotational line total opacity and intrinsic 
FWHM width of individual hyperfine transitions (as fitted by 
the CLASS HFS routine, which assumes equal excitation tem- 
perature for all sublevels), and relative velocity centroids and 
intensities of the main hyperfine groups in our spectra (derived 
from gaussian fits). Note that the hyperfine groups are always 
broader than the intrinsic linewidth derived by CLASS, due to 
optical depth effects and/or to the presence of adjacent hyper- 
fine components too close to be spectroscopically separated. 

To compare with models, beam efficiency correction is an 
important issue. Since L183 (as well as other PSCs) is very 
extended compared to the primary beam of 10"-30", the T m b 
intensity scale is inadequate, as it will overcorrect for source- 
antenna coupling and thus overestimate the true source surface 
brightnes^l (cf. PPABC). To avoid this problem, observed spec- 
tra are corrected only for moon efficiency (~ T^ scale), while 
models are convolved by the full beam pattern of the 30m tele- 
scope (cf. Greve et al.Q]5980, yielding intensities on the same 
T^ scale. For the CSO, we also correct the data for moon effi- 
ciency (77 m0 on ~ 0.8 at 279 GHz) and as the details of the beam 
pattern are not known, the model output is convolved with a 
simple gaussian beam. Because the CSO main beam coupling 
is very good at 1 . 1 mm (t]mb ^ 0.7), the uncertainty of the cor- 
rection should be limited. 

In order to treat correctly the beam coupling to the source, 
we also take into account the fact that the LI 83 core is located 
within a larger-scale N2H + filament oriented roughly north- 
south, as revealed by previous low resolution maps (PPABC). 
This elongation is mostly constant in intensity over ^ 6 arcmin- 
utes. Therefore we approximate the brightness distribution as a 



2 http://www.iram.fr/IRAMFR/GILDAS/ 

3 Teyssier et al. (20021 provide appropriate corrections for 30-m 
data, but only for uniform circular disk sources, and for data taken 
before the 1998 surface readjustment. 

4 The last surface readjustment of the 30-m in 1998 has not been 
characterized in detail, therefore we have scaled down the error beam 
coupling coefficients given in Greve et al. (1998 ) to retrieve the new, 
improved beam efficiencies 



cylinder of length 6', with its axis in the plane of the sky. The 
intensity distribution along the east-west section of the cylin- 
der is taken as the emergent signal along the equator of our 
spherical Monte-Carlo model. We then replicate these values 
along the (north-south) cylinder axis before convolving by the 
antenna beam. For N2D + , the (unpublished) large-scale map 
shows a smaller north-south extent (^ 2') hence we use a cylin- 
der length of 2' for convolving our model. 

In Figs. [1] and [2] we plot the observations in T^ compared 
to our best model (see Section 13.4) . convolved by the antenna 
beam as explained above. 

3. Modelling and discussion 

3.1. Radiative transfer code 

Our spherical Monte-Carlo code is derived from Bernes' code 
(Bernes 1979) and has been extensively tested on models of CO 
emission from dense cores (e.g. Pagani |T998l l. It was recently 
updated to take into account overlap between hyperfine tran- 
sitions occuring at close or identical frequencies. This is done 
by treating simultaneously the photon transfer for all hyper- 
fine transitions of the same rotational line (see also Appendix 
in Gonzalez-Alfonso & Cernicharo 1993). Instead of choos- 
ing randomly the frequency of the emitted photon, a binned 
frequency vector is defined for each rotational transition inside 
which all hyperfine transitions are positioned. All bins are filled 
with their share of spontaneously emitted photons (or 2.7 K 
continuum background photons) and, during photon propaga- 
tion, absorption is calculated for each bin by summing all the 
hyperfine transition opacities at that frequency 



r(v) = J K(v)ds = ^t,0,(v) = J ^icifriy) 



ds 



where k-, is the absorption coefficient, r, the opacity and 0;(v) 
the local frequency profile of the \ ,h hyperfine transition. The 
total number of absorbed photons I (v)e~ T(v) is then redis- 
tributed among hyperfine levels according to their relative ab- 
sorption coefficients 

d/;(v) = — — - X I (v)e Tm 
k{v) 

Overlap is treated for all rotational transitions. Statistical 
equilibrium, including collisions, is then solved separately for 
each hyperfine energy level. 

Collisional coefficients are available for transitions between 
all individual hyperfine levels, but were calculated with He as 
a collisioner instead of H2 (Daniel et al. 120051) . We scaled them 
up by 1.37 to correct for the difference in reduced mass, but 
note that this correction is only approximate (see Willey et al. 
2002 and Daniel et al. 2006a ). The code works for both isotopo- 
logues. As one can probably neglect the variation of the electric 
dipole moment (Gerin et al. 2001}, the N2D + Einstein sponta- 
neous emission coefficients are derived from those of N2H + by 
simply scaling them down by 1.2 3 . Deexcitation collisional co- 
efficients are kept the same as for N2H + (following Daniel et 
al. 12006b) . 
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Fig. 2. N2D + observations (grey 
lines) compared to our best radia- 
tive transfer model for the largest 
N2D + abundance possible (dotted 
histogram in Fig |4^). All data are 
from IRAM 30-m. Only the central 
triplet is represented for the (J: 1-0) 
transition 



The local line profile in our ID code takes into account both 
thermal and turbulent broadening, as well as any radial and ro- 
tational velocity fields. In principle, rotation introduces a devi- 
ation from spherical symmetry. However, for small rotational 
velocities (less than half the linewidth, typically), the excita- 
tion conditions do not noticeably change within a given radial 
shell, as shown by a previous 2-D version of this code (Pagani 
& Breart de Boisanger, 1996). Therefore the (much faster) ID 
version remains sufficiently accurate. 

Our Monte-Carlo model shows that excitation temperatures 
among individual hyperfine transitions within a given rotational 
line differ by up to 15 %, hence the usual assumption of a sin- 
gle excitation temperature to estimate the line opacity (e.g. in 
the CLASS HFS routine) is not fully accurate as already no- 
ticed by Caselli et al. (1995). Neglecting line overlap affects 
differentially the excitation temperature of hyperfine compo- 
nents when the opacity is high enough, mostly decreasing it, 
but also increasing it in a few cases. For example, in our best 
model for the LI 83 PSC, the excitation temperatures of indi- 
vidual hyperfine components change by up to 10%. As shown 
in the last column of Fig. Q] a noticeable difference between 
models with and without overlap appears in the (J:3-2) line 
shape. 

3.2. Grid of models 

Our core model has an outer radius of 13333 AU 2', fixed 
by the maximum width =* 4'of detectable NaLF emission in the 
east-west direction across the PSC, as seen in large scale maps 
(PPABC). We assume that N2I-F abundance is zero outside of 
this radius (due to chemical destruction by CO). 

For all models, the microscopic turbulence is set to 
Av tur b(FWHM) = 0.136kms~ 1 . This was found sufficient to re- 
produce the width of individual hyperfine groups, and is com- 
parable to the thermal width contribution. A small rotational 



velocity field of a few tens of m s _1 was further imposed beyond 
3000 AU (see Table |2]i to reproduce the small anti-symmetric 
velocity shift from PSC center observed at distances beyond + 
30" (cf. Section[2]). The radial velocity was kept to in all lay- 
ers. As seen in Figs. 1 and 2, this approximation already gives 
a remarkable fit to observed line profiles, and is therefore suffi- 
cient to derive the overall temperature and abundance structure. 

We considered a variety of density profiles: single power 
law profiles p oc r _1 , p oc r _1 5 , p oc r~ 2 , as well as broken power 
laws with p oc r _1 out to 4000 AU followed by r~ 2 , or by r~ 3 5 . 
For good accuracy, the density profile was sampled in 330 AU 
= 3" thick shells, i.e. 3 to 10 times smaller than our observa- 
tional beam sizes. In all cases, divergence at r — was avoided 
by adopting inside r = 990 AU the density slope derived from 
the ISOCAM data (Pagani et al. 2004). The density profiles 
were then scaled to give the same total column density N(H2) 
= 1.4 x 10 23 cirT 2 towards the PSC center, as estimated from 
the MAMBO emission map (Pagani et al. 2004). The effect of 
changing the total N(Fb.) on the fit results will be discussed in 
Section I3AT1 

We then considered a variety of temperatures spanning the 
range 6-10 K in the inner core regions. We also tested non- 
constant temperature laws, increasing up to 12 K in the outer- 
most layers and/or in the innermost layers. 

For each combination of density and temperature laws, a 
range of N2FF abundance profiles was investigated. To avoid 
a prohibitive number of cases, abundances were sampled in 6 
concentric layers only, corresponding to the spatial sampling of 
our spectra (except the last, wider layer which encompasses the 
outermost two offsets), and multiple maxima/minima were for- 
bidden. This lead to a total of about 40,000 calculated models 
with 8 free parameters (6 abundances, 1 temperature profile, 
1 density profile). Table [2] lists the values of H2 densities for 
4 density laws, together with the temperature and abundance 
profiles that best fitted the data in each case (cf. next section). 
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Table 2. The 4 main density laws investigated in our models, and the corresponding temperature and abundance profiles giving 
the best fit to N2H + data in each case. The leftmost model gave the best overall fit (see Table [3]) and was used to compute the 
synthetic spectra in Figs. 1 and 2. Its best N2D+ abundance profile is also listed. Parameters are held constant inside a given 
radial shell, except for the rotational speed which is linearly interpolated. 2.34(6) reads 2.34 x 10 6 . 







P 


oc r 1 


r~ 2 (Best model)" 




p oc 


r 1 




p oc r 


-1.5 




p oc 


r- 2 


Radius 


V. rot. 


Density 


T 

1 kin 


Abundance 


Abundance 


Density 


T 

1 kin 


Abundance 


Density 


T 

1 kin 


Abundance 


Density 


T 

1 kin 


Abundance 


A TT 


(Kill S ) 


(cm ) 




N 2 H+ 


N 2 D+ 


(.cm ) 




NT T 1 1 
IN 2 1 1 


(cm ) 


(A) 


m rr+ 
lN 2 rl 


(cm ) 




M T 1 ! 
,N 2 1 1 


330 





2.34(6) b 


7 


2.4(-ll) 


2(-H) 


2.34(6) b 


7 


7(-ll) 


3.30(6) c 


7 


l(-ll) 


3.46(6) c 


9 


l(-ll) 


660 





2.05(6) b 


7 


2.4(-ll) 


2(-H) 


2.05(6) b 


7 


7(-H) 


2.89(6) c 


7 


K-ll) 


3.03(6) c 


9 


l(-ll) 


990 





1.55(6) b 


7 


8.5(-ll) 


4(-H) 


1.55(6) b 


7 


7(-H) 


2.19(6) c 


7 


7(-H) 


2.29(6) c 


9 


3(-ll) 


1320 





1.16(6) 


7 


8.5(-ll) 


4(-ll) 


1.16(6) 


7 


7(-ll) 


1.19(6) 


7 


7(-ll) 


1.46(6) 


9 


3(-H) 


1650 





9.27(5) 


7 


8.5(-ll) 


4(-H) 


9.27(5) 


7 


7(-H) 


7.74(5) 


7 


7(-H) 


9.09(5) 


9 


3(-ll) 


1980 





7.73(5) 


7 


8.5(-ll) 


4(-ll) 


7.73(5) 


7 


7(-H) 


5.55(5) 


7 


7(-H) 


6.10(5) 


9 


3(-H) 


2310 





6.62(5) 


7 


1.1(-10) 


3(-H) 


6.62(5) 


7 


1(-10) 


4.20(5) 


7 


4(-10) 


4.39(5) 


9 


K-10) 


2640 





5.80(5) 


7 


1 If— 10) 


3(-ll) 


5.80(5) 


7 


1(-10) 


3.34(5) 


7 


4(-10) 


3.54(5) 


9 


1(-10) 


2970 





5.15(5) 


7 


1.1(-10) 


3C-11) 


5.15(5) 


7 


1(-10) 


2.74(5) 


7 


4(-10) 


2.68(5) 


9 


K-10) 


3300 


0.01 


4.64(5) 


7 


l.l(-lO) 


3(-H) 


4.64(5) 


7 


1(-10) 


2.29(5) 


7 


4(-10) 


2.20(5) 


9 


K-10) 


3630 


0.01 


4.21(5) 


7 


1.5(-10) 


3(-H) 


4.21(5) 


7 


1.5(-10) 


1.95(5) 


7 


3(-10) 


1.83(5) 


9 


6(-10) 


3960 


0.01 


3.54(5) 


7 


1.5(-10) 


3(-H) 


3.86(5) 


7 


1.5(-10) 


1.70(5) 


7 


3(-10) 


1.46(5) 


9 


6(-10) 


4290 


0.02 


3.01(5) 


7 


1.5(-10) 


3(-H) 


3.57(5) 


7 


1.5(-10) 


1.49(5) 


7 


3(-10) 


1.22(5) 


9 


6(-10) 


4620 


0.02 


2.60(5) 


7 


1.5C-10) 


3(-H) 


3.31(5) 


7 


1.5(-10) 


1.32(5) 


7 


3(-10) 


1.07(5) 


9 


6(-10) 


4950 


0.02 


2.26(5) 


7 


1.3C-10) 


5(-12) 


3.09(5) 


7 


5(-H) 


1.18(5) 


7 


3(-10) 


9.27(4) 


9 


4(-10) 


5280 


0.05 


1.99(5) 


7 


1.3(-10) 


5(-12) 


2.90(5) 


7 


5(-H) 


1.07(5) 


7 


3(-10) 


8.5(4) 


9 


4(-10) 


5610 


0.05 


1.76(5) 


7 


1.3(-10) 


5(-12) 


2.73(5) 


7 


5(-H) 


9.67(4) 


7 


3(-10) 


7.20(4) 


9 


4(-10) 


5940 


0.05 


1.57(5) 


8 


1.3(-10) 


5(-12) 


2.58(5) 


8 


5(-ll) 


8.82(4) 


7 


3(-10) 


6.34(4) 


9 


4(-10) 


6667 


0.05 


1.27(5) 


8 


1(-10) 


< 4(-12) 


2.02(5) 


8 


2.5(-ll) 


8.10(4) 


7 


1.5(-10) 


5.25(4) 


9 


4(-10) 


7733 


0.05 


8.50(4) 


9 


1(-10) 


< 4(-12) 


1.57(5) 


9 


2.5(-ll) 


6.77(4) 


7 


1.5(-10) 


4.55(4) 


9 


4(-10) 


8867 


0.05 


6.50(4) 


10 


1(-10) 


< 4(-12) 


8.50(4) d 


10 


2.5(-ll) 


5.42(4) 


7 


1.5(-10) 


3.46(4) 


9 


4(-10) 


10000 


0.02 


5.20(4) 


11 


1(-10) 


< 4(-12) 


5.20(4) d 


11 


2.5(-ll) 


4.41(4) 


7 


1.5(-10) 


2.68(4) 


9 


4(-10) 


13333 





3.45(4) 


12 


1(-10) 


< 4(-12) 


3.45(4) d 


12 


2.5(-ll) 


3.68(4) 


7 


1.5(-10) 


1.70(4) 


9 


4(-10) 



a p oc r 1 (r < 4000 AU), p oc r 2 (r > 4000 AU) 

b The first three layer densities do not follow a power law but the density profile derived from the ISOCAM absorption map (Pagani et al. 
|2004t 

c Same as note (b) above, but rescaled to maintain a constant total column density of 1.4x 10 23 cirT 2 (see text). 

d For the pure p oc r _1 case, density is not falling off fast enough to reach ambient cloud density and we correct for this in the last 3 layers. 



3.3. Selection criteria 

Our best fit selection criteria are based on the use of several 
X 2 and (reduced) xl estimates. Because the N2H + (J: 1-0) data 
have much more overall weight than the (J: 3-2) data (more iso- 
lated hyperfine groups per spectrum, more observed positions, 
higher signal-to-noise ratio), computing a single xl P er model 
is hardly sensitive to the quality of the (J:3-2) line fit. However, 
the (J: 1-0) spectra being optically thick, they give little con- 
straints on the temperature and N2H + abundance of the inner- 
most layers (within a 6" radius). Those are only constrained by 
the IRAM (J:3-2) observations. The CSO (J:3-2) spectrum at 
(24", 0") also adds interesting information on the density pro- 
file. 

To constrain the models we thus evaluate the goodness of fit 
independently for the (J: 1-0) and (J:3-2) data, by computing ;f 2 
values on 4 types of measurements : 1) full area of each of the 7 
(J: 1-0) hyperfine components, at each of the 7 offset positions 
(49 values). This measurement is sensitive to the global temper- 
ature and abundance profiles. Fitting separately each hyperfine 



component allows to discard solutions with overall identical 
area but different hyperfine line ratios; 2) area of the 5 central 
channels (= 0.15 kms 1 ) of each hyperfine component for the 
same lines (49 values), in order to reject self-absorbed profiles, 
i.e. line profiles which have the same total area but are wider 
and self-absorbed; 3) total area of the (J:3-2) CSO spectrum at 
(24", 0") (1 value). This single measurement helps to constrain 
the density profile; 4) total area of the (J:3-2) 30-m spectrum 
(1 value) which constrains the abundance (and the temperature 
to some extent) at the core center. These four measurements 
are of course not completely independent from each other, and 
trying to improve one fit by changing a model parameter can 
degrade one or several of the others. We used xl f° r the first 
two (with 41 degrees of freedom) and^ 2 for the last two (only 
one measurement) so that in all cases Ay 2 « 1 is equivalent to 
a lcr deviation. With this choice, a 3cr variation is equivalent to 
Axl ~ 1.5 for the first two and Ay 2 = 9 for the last two. 

After having run all of the models in our grid, we first used 
a single, global xl combining the four types of measurements 
to visualize the global fit quality of the various models. Fig. [3] 
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Fig. 3. Global;^;; for all density profiles and temperature cases. 
An arrow indicates that the majority of the model layers is at 
the given temperature but that outer layers have rising temper- 
atures. Only one (circled) solution satisfies all 4 individual^ 2 
fitting (see text) 



plots the global^- 2 as a function of temperature, for the various 
density laws that we investigated. Symbols with arrows indi- 
cate models where temperature rises in the outermost layers. 
For each density-temperature pair, we plot only the smallest ;f 2 
obtained by adjusting the abundance profile. 

It is readily seen that the best fits are found for relatively 
low overall core temperatures ^ 6.5-8 K. Hence, the (J: 1-0) 
data alone (which dominate the global xl) set a strong con- 
straint on the core temperature: indeed, the combination of 
strong opacity and relatively low intensity in this line requires 
low excitation temperature. Because the (J: 1-0) line should be 
thermalized in the dense center, this rules out kinetic tempera- 
tures above 8 K (for our adopted total column density of 1 .4x 
10 23 cm" 2 ). 

To further discriminate among models, we then looked in- 
dividually at the four quality indicators described above, in par- 
ticular those related to the (J:3-2) lines. Table[3]lists these 4 val- 
ues as well as the global xl f° r various models selected from 
Fig. [3] namely: the best model fits for a given constant core 
temperature (from 6 to 10 K), and the best model fits for a 
given density law (with corresponding temperature and abun- 
dance distributions given in Table[2]). 

It can be seen in Table [3] that for Tki n = 10 K the (J:3- 
2) lines are not well fitted (4 and 5 cr deviations for the CSO 
and IRAM lines respectively), hence such a high temperature 
seems ruled out here. Overall, it is rather difficult to have low 
X 2 values in both the (J: 1-0) and the (J: 3-2) indicators. Only 
one density-temperature combination (indicated by bold faces 
in Table [3) reaches low values in all 4 individual xl an d X 2 > 
each within 1 cr of their minimum value, among all the models 
in Table[3] This combination is referred to as our "best model" 
in the following, and its inferred physical parameters are dis- 
cussed in the next section. It can also be noticed in Table[3]that 
the model with T^n = 7 K has a slightly better global xl than 
the best model (2.1 compared to 2.3) but is to be eliminated as 
the (J: 3-2) lines are 2 and 3.5 cr off the mark. 



3.4. Physical conditions in the L183 main core 

3.4.1 . Density and temperature structure 

Our best density-temperature combination (circled 5-branch 
star symbol in Fig. [3j has a density law of r _1 out to r « 4000 
A.U., consistent with the ISOCAM profile (Pagani et al. 2004) 
and a slope of r~ 2 beyond. The temperature is constant at 7 K 
out to 5600 A.U. and increases up to 12 K at the core edge 
(see Table Further increasing the number of warmer lay- 
ers failed, as well as introducing warmer layers in the center. 
Keeping the temperature constant at 7 K everywhere signifi- 
cantly worsens the fit of the (J:3-2) lines (see TableO. 

The main uncertainty in our temperature determination 
stems from the assumed total gas column density through the 
core. Decreasing/increasing it by a factor 1 .4 (the typical un- 
certainty found by Pagani et al. 2004) changes all volume den- 
sities by the same factor. To recover the same N2H + emission, 
we find that the kinetic temperature in our models must be in- 
creased/decreased by ^1 K (respectively, while the N2H + abun- 
dance must be scaled accordingly to keep the same column 
density). Since the dust temperature in the L183 core is 7.5 + 
0.5 K (Pagani et al. 2004), gas and dust are thermalized inside 
this core within the uncertainties. 

A different result was found by Bergin et al. (12006 ) in the 
B68 PSC, where outer layers emitting in CO are consistent with 
a kinetic temperature of 7-8 K while NH3 measurements indi- 
cate a higher temperature of 10-1 IK in the inner 40". This 
inward increase in gas temperature was attributed to the lack of 
efficient CO cooling in the depleted core center, and requires 
an order of magnitude reduction in the gas to dust coupling, 
possibly due to grain coagulation (Bergin et al. 2006). Such an 
effect is not seen in LI 83, despite strong CO depletion within 
the inner V- 6600 AU radius. Our best fit temperature law is 
more consistent with the thermo-chemical evolution of slowly 
contracting prestellar cores with standard gas-dust coupling co- 
efficients, which predicts low gas temperatures ^ 6 K near the 
core center, increasing to ^ 14 K near the core surface (Lesaffre 
etal.|2005]l. 

3.4.2. N 2 H + abundance profile 

With our best density-temperature model, several N 2 H + abun- 
dance profiles give equally good fits (within <lcr on all 4 qual- 
ity indicators). From all these acceptable abundance profiles, 
we determined a median profile (plotted as a solid histogram 
in Fig. UK), with error bars representing the range of acceptable 
values in each layer (though not any combination of these val- 
ues does fit the observations). This median profile, whose val- 
ues are listed in Table|2] is used to compute the fit displayed in 
Fig.Q] The total N 2 H + column density is 1.2+0.1 x 10 13 cm -2 , 
comparable with the value reported by Dickens et al. (|2000) 
and a factor of ~2 below that in Crapsi et al. d20051 >. 

The maximum N 2 H + abundance is 1.5*qj x 10~ 10 , the 
same as found by Tafalla et al. d2004t in L1498 and L1517. 
However, there is a definite drop in abundance in the inner lay- 
ers. This drop is imposed by the fit to the (J:3-2) IRAM central 
spectrum (9" beam), which is quite sensitive to variations in 
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Table 3. Quality of fit evaluated individually on four different sets of measurements, and global for selected models from 
Fig.[3](see text). A difference of 1 corresponds to lcr in all cases. Values for the best model are in bold face. The last rows show 
the effect of changing only the central N2H + abundance in this best model 



(J: 1-0) (J:3-2) 
model parameters Total area xl 5 channels area xl CSO Total area x 1 IRAM Total area x 1 Global^ 

Best fit for each temperature (adjusted in abundance) 



T = 6K (par',r 2 ) 4.7 3.8 1.5 1.0 3.9 

T = 7K (pccr-V 2 ) 1.9 2.4 12.6 4.1 2.1 

T = 8K (peer- 15 ) 3.1 3.2 0.2 3.4 2.9 

T = 9K (peer 2 ) 4.5 5.6 7.5 6.8 4.7 

T=10K (peer 2 ) 4.8 6.3 16.9 28.7 5.5 



Best fit for each density law (adjusted in abundance) 

peer" 1 ,!-- 2 (T = 7^12K) 2.4 2.7 0.3 0.0 2.3 

peer 1 (T = 7->12K) 3.7 3.8 1.2 2.9 3.4 

peer 15 (T = 7K) 3.4 3.2 0.1 0.1 2.9 

pocr- 2 (T = 9K) 4.5 5.6 7.5 6.8 4.7 



Varying the central N 2 H + abundance at r < 660 AU in the best model 

X(N 2 H + )= 10~ 12 2.6 2.8 0.3 2.0 2.5 

X(N 2 H + )= 8xl0- 12 2.6 2.8 0.3 1.0 2.5 

X(N 2 H)= 2.4x10 11 2.4 2.7 0.3 0.0 2.3 

X(N 2 H + )= 8x10-" 2.1 2.4 0.5 9.2 2.1 



radius (arcsec) radius (arcsec) 

10 1 10 2 10 1 10 2 




radius (AU) radius (AU) 



Fig. 4. a N2D + and N2H + abundances for the best model. For N2H + , the median abundance profile is plotted (the corresponding 
model output is displayed in Fig. 1). Error bars indicate the range of possible values (see text). For N2D + , a range of values is 
given: The upper, dotted histogram gives the best fit to the (J:2-l) and (J:3-2) lines, and was used to compute the model in Fig. 2. 
The dashed histogram gives the best fit to the (J: 1-0) and (J:2-l) lines, b density profile of the best model and N2D + /N2H + ratio 
range 



the central abundance of N2H + , as illustrated in the last rows 
of Table [3] In contrast, the (J: 1-0) line is optically thick at the 
core center and the^J is thus not sensitive to this parameter (cf. 
Tabled. 

Exploring a broad range of abundance profiles, we find only 
a limited range of possible N2H + abundances in the inner layer 
of 6" radius (r < 660 AU): values above 4.5 x 10~ n always 
produce too strong (J:3-2) emission in the central IRAM beam 
(X 1 > 1), while abundances below 10~ 12 give a (J:3-2) line 
marginally too weak. We thus derive a range of 2.4 x 10~ n 
for the central N2H + abundance. However, we note that our 



results for N2D + would favor a value of at least 10 11 (see next 
section). 

Compared with the maximum abundance, this gives a vol- 
ume depletion factor of 6^ 3 at the core center. The median 
abundance drop is less than (but marginally consistent with) 
that expected from simple geometrical arguments in PPABC, 
but it confirms that the leveling of N2H + intensity seen across 
the dust peak is not due to pure opacity effects. As can be seen 
from Fig. Hk-b, depletion starts at a density of 5-7 x 10 s cm -3 
and increases as density goes up, in agreement with the con- 
clusions of PPABC. The abundance also drops slightly in the 
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outermost regions, possibly due to partial destruction by CO. 
Indeed, the outermost layers of the model have densities of a 
few 10 4 cm -3 which is the limit above which CO starts to de- 
plete in this source (PPABC). 



3.4.3. N 2 D + abundance profile 

As seen in FigfJ] good fits to the N2D + data may be obtained 
with the same density and temperature profile as our best model 
for N2H + , although it is difficult to reproduce simultaneously 
the intensities of both (J: 1-0) and (J:3-2) lines. We plot in Fig. 
[4j the N2D + abundance profiles giving the best fit either for 
(J: 1-0) and (J:2-l) (dashed histogram) or for (J:2-l) and (J:3-2) 
(dotted histogram). The latter was used to produce Fig. [2] A 
better fit could be obtained with a temperature of 8 K instead 
of 7 K. This might indicate that the collisional coefficients with 
He are somewhat inaccurate to represent those with H2. Indeed, 
it has been shown for NH3 that collisional coefficients with He 
could differ by a factor up to 4 with respect to those computed 
with para-H2 (Willey et al. 120021) . A similar problem may occur 
here (see also the discussion in Daniel et al. l2006al) . Therefore 
we consider that the temperatures are compatible within the 
uncertainties on the collisional coefficients. 

The N2D + abundance profile is quite different from that 
of N2H + (cf. Fig. Hk). Its abundance is essentially an upper 
limit beyond 6000 A.U. It increases sharply by about an order 
of magnitude in the region between 600 and 4000 A.U., then 
slightly drops by a factor of 2-2.5 in the core center, reach- 
ing an abundance between 1.3 and 2x 10~ n . The low optical 
depth of the line (t = 0.84 for the strongest hyperfine com- 
ponent, J FF ' : I23— O12) allows to measure the contribution of all 
layers to the emission and to determine with relatively little 
uncertainty the abundance profile. In particular, for the chosen 
density and temperature profiles, we find that the abundance of 
N2D + in the center of the cloud cannot be below 10 _1 1 . Hence, 
we can set tighter constraints on the N2D + abundance at the 
core center than was possible for N2H + . 

The N2D + /N2H + ratio obtained by comparing the N2D + 
range of abundances to the median N2H + abundance profile 
is plotted in Fig. @J3. The deuteration ratio varies from an up- 
per limit <0.05-0. 1 away from the core to a very high factor of 
0.7±0.12 in the depletion region. This is larger than what Tine 
et al. (2000) reported 48" further north in the same source, and 
3-4 times larger than what Crapsi et al. {2005) report towards 
the PSC. Deuterium enrichment is thus very strong and can cer- 
tainly be linked to the strong and extended H2D + line detected 
towards this source (Vastel et al. 2006). If we consider the full 
range of possible values for N2H + itself, the deuteration ratio 
in the inner layer ranges from ~0.3 to > 20. However, as D2H + 
has not been detected in this source despite its strong H2D + line 
(Vastel et al. 2006), an enrichment above 1 seems unprobable, 
suggesting that the central abundance of N2H + is probably at 
least lCT 11 . 
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Fig. 5. GBT NH3 (1,1) inversion line towards the reference po- 
sition with CLASS NH3 fit (shifted by +3 K for clarity) and 
fit residual (shifted by -0.6 K). The fit yields a total opacity of 
24.2 (±0.4), and a global velocity of 2.3672 (±0.0002) km s" 1 



3.5. Comparison with NH3 temperature estimates 

The low kinetic temperatures <8 K inferred from our N2H + and 
N2D + modelling are somewhat lower than previous tempera- 
ture estimates in L183 from NH3 inversion lines, which gave 
values in the range 9-10 K (Ungerechts et al. [T980, Dickens et 
al. 120001 up to 12 K (Swade et al. [19891 . However, these NH 3 
spectra were not obtained towards the PSC center itself, and 
had relatively low angular resolution for the last two. We thus 
briefly reconsider this issue using our NH3 spectra obtained 
at the PSC center position, which also benefit from the higher 
resolution and beam coupling of the new GBT. The NH3 spec- 
tra were analyzed in the standard way, as discussed by Ho & 
Townes ( 1983 ) and Walmsley & Ungerechts (119831) . Using the 
CLASS NH3 fitting procedure, we found a total opacity of 24 
for the (1,1) inversion line, and an excitation temperature of 5.5 
K assuming a beam filling-factor of 1 (Figs.|5]&|6]l. Making the 
usual assumption of constant temperature on the line of sight 
and of negligible population in the non-metastable levels, the 
intensity ratio of the (2,2) to (1,1) main lines indicates a ro- 
tation temperature of 8.4 ±0.3 K which should correspond to 
a kinetic temperature of 8.6 ±0.3 K. Hence, NH3 emission to- 
wards the PSC center indicates an only slightly higher kinetic 
temperature than N 2 H + (8.6 ±0.3 K instead of 7 ±1 K), almost 
equal to that obtained with N2D + (~8 K). The discrepancy be- 
tween NH3 and N2H + is thus not as large as originally thought. 
Both tracers point to very cold gas in the core, close to thermal 
equilibrium with the dust. 

Reasons for a possible difference between NH3 and N2H + 
temperature determinations in LI 83 include the following: 1) 
the higher sensitivity of NH3 to the warmer outer layers of the 
core, since NH3 inversion lines are much easier to thermalize 
(n cr j t ~ 2000 ctrT 3 ) than N2H + lines. 2) a slight overestimate 
in the total column density towards the PSC (reducing it to 
10 23 crrT 2 , the temperature has to be raised to 8 K to compen- 
sate for the density decrease and recover the same N2H + emis- 
sion). 3) systematic errors introduced by the use of collisional 
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Fig. 6. GBT NH3 (2,2) inversion line towards the reference po- 
sition. A gaussian fit to the main line indicates a velocity of 
2.362 (+0.005) km s~' 

coefficients with He instead of H2 forN2H + (see Sect. l3.4.3l . 4) 
concerning NH3, standard hypotheses leads to a puzzling dis- 
crepancy between T ex = 5.5 K and T rot = 8.4 K towards the 
L183 PSC center, where NH3 inversion lines should be fully 
thermalized. Beam dilution has been invoked for giant molec- 
ular clouds (thus increasing T ex ), but it seems unreasonable to 
extend this to dense cloud cores (Swade 1 1 9891 and references 
therein). A Monte-Carlo code (or equivalent) would be needed 
to model NH3 taking into account the strong density gradients 
present in PSCs, and the possible population of non-metastable 
levels at very high densities. 

4. Conclusions 

1. We have presented a new Monte-Carlo code (available 
upon request to the author) to compute more realistically 
the NLTE emission of N2H + and N2D + , taking into account 
both line overlap and hyperfine structure. This code may be 
used to infer valuable information on physical conditions in 
PSCs. 

2. The best kinetic temperature to explain N2H + observations 
of the LI 83 main core is 7+1 K (and ~ 8 K forN 2 D + ) inside 
5600 AU, therefore gas appears thermalized with dust in 
this source. 

3. There is no major discrepancy with NH3 measurements 
which also indicate very cold gas (8.6 ±0.3 K) towards the 
PSC. 

4. We have found a noticeable depletion of N2H + by a fac- 
tor of 6^ 3 , and of N2D + by a smaller factor of 2-2.5. This 
smaller depletion is probably due to a strong (0.7+0.12) 
deuterium fractionation, consistent with the detection of 
H2D + in this core. 

5. N2D + should be a useful probe of the innermost core re- 
gions, thanks to its low optical depth combined with its 
strong enhancement. 
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